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We derive new analytical results for the hydrodynamic force exerted on a sinusoidally 
oscillating porous shell and a sphere of uniform density in the Stokes limit. The coupling 
between the spherical particle and the solvent is done using the Debye-Bueche-Brinkman 
model, i.e. by a frictional force proportional to the local velocity difference between 
the permeable particle and the solvent. We compare our analytical results and existing 
dynamic theories to Lattice-Boltzmann simulations of full Navicr-Stokcs equations for 
the oscillating porous particle. We find our analytical results to agree with simulations 
over a broad range of porosities and frequencies. 
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1. Introduction 

The degree of porosity affects sedimentation and aggregation dynamics of suspended 
particles. For instance, high-porosity, low-density particles have been found suitable for 
efficient delivery of therapeutics into the systemic circulation through inhalation. Such 
treatment is made possible by careful engineering of porous particle structur e and dy- 
nami cs to circumvent pulmonary mechanisms for removing deposited particles ( Edward's] 
19971) . Besides clinical applicability, models of fluid flow through porous media have 



been developed and tested, for instance, to improve t he efficiency of oil recovery by 
fluid injection (jBabadaglil 120031) and ultrasonic waves (lAmro &; Al- Homadhi 2006), to 



characterize structural properties in pulp and paper science ( Ramaswamv et al. 2004) 



and to identify conditions that cause colloid detachment from surfaces in porous me- 
dia ( Bergendahl k, Grassollioooh . The response to the surrounding flow depends on the 
mass distribution within the porous particle. Steady-state response to hydrodynamic 
forces and torques is well understood, but the dynamics of permeable particles still 
poses several unanswered questions. Only recently has it become p ossible to study diffu - 
sive properties of concentrated suspensions of permeable particles (jAbade et a/.ll2010al ). 

Porosity together with understanding of hydrodynamic forces in a corrug ated nanochan- 

nel co uld be used to ta ke advantage of size-dependent transport properties ( Del Bonis- O'Donnell et al 
|2009( ) and separation ( Fu et al. 2007 ) of nanofluids ( Sparreboom et al. 2010l ). 

The first a nalytical results o n stead y-sta te drag on a perme able homogeneous sphere 
date back to iDebve fc Buechd (|l948t) and iBrinkman l dl947allbh . Since then, numerous 
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work s (lFelderholl975l:lFelderhof fc Deutchlll975l : lbeutch fc Felderhodl975UBhatt fc Sacheti 
19941 : ICichocki fc Felderhofl l2009h have been published on drag force and torque on 
spherical particles of d ifferent mass distributions and internal structure. Fo r example, 



Bhat t fc Sachetil (| 1994T ) investigated a porous shell of finite thickness. Recently. ICichocki fc Felderhof 
( 20091 ) solved a related problem where the shell was wrapped around a solid core. So far, 
studies on the subject have largely comprised theoretical calculations in the Stokes ap- 
proximation of the Navier-Stokes equation. 

In this work, we first compare the results of such calculations of steady-state quantities 
to computer simulations of the Navier-Stokes equations by the well-established Lattice- 
Boltzmann method (LB). We will show that our simulations give quantitative agreement 
with theoretical predictions without any adjustable fitting parameters at all levels of 
permeability for the steady-state case. 

We then examine the dynamic case in which addit ional complicatio n s aris e as the 
particle moves in the fluid in an oscillatory manner. Looker fc Carniei ( 2004 ) used a 



perturbative expansion to find the hydrodynamic force on a slightly permeable sphere. 
They found significant differences in the fluid velocity aro und and in the hydrodynamic 
force on the particle in a frequency ra nge from 1 to lOMHz. lVainshtein fc Shapiro! (|2009h 
generalized the original Stokes ( 190lh result of the hydrodynamic force on a sinusoidally 
oscillating solid particle with a no-slip boundary condition. They formulated the problem 
such that changes in both the velocity and acceleration-dependent part of the dynamic 
force could be quantified as the frequency of the oscillation, the porosity of the particle or 
as the boundary condition on the surface of the particle was changed. In the high-porosity 
limit, where the Brinkmann (3 parameter tends to zero, th e hydrodynamic drag force on 
the particle should approach zero. However, the model of Vainshtein fc Shapiro! ( 20091 ) 
gives a finite hydrodynamic drag force on the particle for all values of f3. In this work, we 
present a new analytical derivation for both a porous shell and a uniform density porous 
sphere. Our results give a physically consistent hydrodynamic force for all values of the 
coupling parameter. We then find our new result to provide better overall agreement with 
simulations of a particle oscillating in the fluid. Our tests are performed at 0.06-28 MHz 
for particles of radii between 80 and 700 nm. 



2. Model 

2.1. Time- dependent Stokes, Darcy and Brinkman equations 

We use the Dcbye-Buechc-Brinkman (DBB) model, which allows one to study generic 
hydrodynamic effects between a solvent and porous particles with few parameters. Typi- 
cally, this model has been studied theoretically in the steady-state case. However, we are 
interested in oscillating particles so will look at the time-dependent model. The coupled 
porous particle-fluid system in the case of an incompressible fluid in the small Reynolds 
number (Re) limit, can be described by the time-dependent linearized Navier-Stokes 
equations 

V-u = 0; (2.1a) 
p<9 t u = ? ? V 2 u- Vp + f, (2.16) 

where p is the fluid mass density, u the fluid velocity, rj the shear viscosity and p is the 
pressure. Here, the presence of the porous particle is characterized by the force density, 

f = 7n(r)(v-u); (2.2a) 
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where the coupling constant 7 has units mass per time, v is the local velocity of the 
particle at the point r, which contains contributions from centre of mass and rotational 
motion. The "node" density n(r), which has units of inverse volume, has a constant value 
A inside the particle and zero outside the volume B(t) of the particle. Outside B(t), where 
f = 0, (|2.ip is commonly referred to as the unsteady Stokes equation (it is missing the 
nonlinear pu ■ Vu term present in the full Navier-Stokes equations). 

Inside B(t), t he fluid flow inter ac ts directly with the nod es and (|2.ip is referred to as the 
DBB equation (|Brinkmanll947allbl ; [Debve fc Buechelll948l) . The shape of the particle B{t) 



can be varied to gi ve a shell, uniform- density sphere or other distribution. The product 
7A is equal to r]n (|Abade et alj|2010al) . where K 1 is the hydrodynamic screening length 



and k~ 2 = 77/ (7A) is the constant permeability of the particle. 

The DBB equation is a mean-field description of fluid flow in the porous particle under 
the assumption that the particle radius R is large enough compared to the mean pore size 
~ y/r)/ (7A). If one further neglects the Laplacian term ?7V 2 u inside B{t), then one 
arrives at the Darcy model, for which the fluid velocity inside the particle is independent 
of r = |r|. The viscosity is assumed to be rj both inside and outside of B(t). 

In the frame of reference with the origin at the centre of mass of our porous particle, we 
define a spherical coordinate system (r, ip, 9) via (x, y, z) = r(sin 9 cos ip, sin 9 sirup, cos 9). 
As we only consider axisymmetric flows, the solution to (|2.ip is independen t of ip, and 
u, p and the fluid stress tensor cr will only depend on (r,6) (jGraebell l2007t ) . Once the 
dependence of a = cr(r,9) is known, one may proceed to calculate the hydrodynamic 
for ce F and torque T acting at the centre of mass of our porous particle from the stress 
as ( Landau fc Lifschitzl[T987h 



F = / cr -e r dS; (2.3) 

JdB(t) 



T = r e r x cr ■ e r dS, (2.4) 

JdB(t) 

where (e r ,e v ,eg) are the unit basis vectors in spherical coordinates and dB(t) is the 
boundary of B(t). 

The porosity-dependent force and torque exerted on the particle by the fluid can also 
be calculated directly using Newton's third law from the integral of the negative of the 
force density on the fluid: 



F = / -fd 3 x= / - 7 n(r)(v - u)d 3 x 

JB(t) JB(t) 

= -/ 77i(r)(v cm + wx(r-r cm )-u(r))d 3 a:. (2.5) 

JB(t) V ' 

where w is the angular velocity of the particle and v cm its centre-of-mass velocity. 
Schematics of the spherical particles we consider arc shown in figure[TJ The hydrodynamic 
torque T on the particle then reads 

T = / (r - r cm ) x (-f ) d 3 x 

JB(t) 

= -/ 7n(r)(r - r cm )( v cm + w x (r - r cm ) - u(r)J . (2.6) 

JB(t) v ' 

In this work, the particle is either held fixed, v cm = 0, or its velocity is set to be 
sinusoidal along the z axis. We separate t he analytical solu tion to (|2.ip to parts inside 
and outside B(t) following conventions of Felderhofl (1975). The parts are matched as 
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Figure 1. Schematic of the (a) uniform-density sphere with node density ns P (r), (b) infinitely 
thin shell nsh(r) and (c) a shell of finite thickness/annulus riAn(r) which is hollow on the inside 
where the Stokes approximation applies. 

detailed in sections 13.31 and 13.41 by requiring the velocity and stress fields to be equal on 
the boundary of B(t). 

2.2. Previous work 

A few other authors have also examined the time-dependent case. lLooker fc Carniel(l2004h 



performed a perturbative expansion to find the force exerted on a rigid, weakly perme- 
able sphere of radius R oscillating in an incompressible fluid. The dimcnsionless pertur- 
bation parameter e = (k_R) _1 was studied in the range [0,0.05], and they consistently 
neglected terms of order C(e 2 ) and higher. They modeled the porous sphere by applying 
the Beavers- Joseph-Saffman (BJS) boundary condition on its surface to order e. They 
solved the homogenized unsteady Stokes equation by assuming the particle to be imper- 
meable enough s o that the flow externa l to B (t) cannot penetrate the particle interface. 
More recently, IVainshtein fc Shapiro! ( 20091) also studied an oscillating sphere of uni- 



form permeability. They used (|2.ip outside the sphere (as f(r) = for r ^ B(t)) and 
inside the sphere, their dynamic equation, 

pd t u = ??V 2 u - Vp - ?7K 2 u, (2.7) 

did not contain the particle velocity. This equation corresponds to (|2.ip and (|2.2[) inside 
B(t) only if v = in (|2.2[) . Instead of introducing the particle-fluid interaction in (|2.7[) , 
they introduced it boundary condition 

u(r = Re r ) - v = u(r = Re r ), (2.8) 

where the particle velocity is v = v^e lut e Z: u refers to the fluid velocity inside B(t) and 
u to that outside B(t). They also required continuity of components of the stress tensor 
o~ rr and a r g at r = Re r when the Laplacian term was included in (|2.7[) and continuity 
of pressure when it was neglected. Note that boundary condition (|2.8[) is not equivalent 
to the DBB model. Any attempt to remove the particle velocity in the equations of 
motion of the DBB model (which enters in the f in (|2.1[0 would end up introducing the 
particle inertia in (|2.1[) (arising from the equation of motion for the particle). Vainshtcin 
and Shapiro's model has neither v nor dw/dt in (|2.T[) . We find it difficult to motivate 
Vainshtein and Shapiro's model physically and will demonstrate that (|2.7p and (|2.8p 
actually give a physically incorrect limit for the time-dependent case for small (3. 

2.3. Simulation method 

We will compare the Stokes theory results to simulations of the full Navier-Stokes equa- 
tions through the well-established LB method. The mass and momentum conservation in 
a flui d are expressed at the Navier-Stokes level as ( Batchelor 1967 ; Landau fc Lifschit3 



1987 ) 



d tP + d a (pu a ) =0 (2.9) 
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and 



d t (pu a ) + dp{pu a up) = -d a p + f a 

+0/3 (rj(d a Up + 8pU a - -drjU^Sa^j + (djUySaP^J , 



(2.10) 



where r\ and £ are the shear and bulk viscosities and p is the fluid pressure. In this work 
we will use a pressure with linear dependence on density, i.e. p = pv^Sap, where v s is the 
speed of sound. This can be viewed as an ideal gas equation of state or the first term in 
a Taylor expans ion of the pressure about fixed density in which case v% is the isentropic 
compressibility ( Kell 1970f ) . The coupling to the particle phase appears through the force 
density f a . 

Our lattice Boltzmann (LB) fluid algorithm, summarized in Appendix lAl repro duces 
dH9l) and (|2~T0|) in the form typical to most LB algorithms (|Chen fc Doolenlll998l ). The 
shear viscosity in the model is y = prvj:. /3, where v c = Ax/ At is a lattice velocity, and 
C = ?](5/3 - Zvl/vl) (jSwift et al.lll995h . In this paper, r will be chosen in all cases so 



that i] = 0.02gcm s , twice the viscosity of water, and p = lgcm . The speed of 
sound v s is chosen to be = (v s < v c is required for stability in LB algorithms). 
This is sufficiently large so that the fluid is approximately incompressible for steady-state 
situations (largest variation in p < 0.1%). In this work, unless stated we use a time step 
At = 1 ns and a mesh resolution of Ax = 100 nm. 
In the simulation, the node density is discrete, 



N 



(2.11) 



That is, the particle is an extended spherical object consisting of N nodes/constituents 
positioned at and moving at velocity Vj around the centre-of-mass coordinate r cm 
whose velocity is v cm . The nodes are coupled to the fluid lattice locally by weighted inter- 
polat ion, which has been used for polymers consisting of p oint particles ( Ahlrichs &: Diinwe 
19981 ). a nanowire immersed in a nematic liquid crystal ( Smith fc DennistorJ 120071 ) and 
most recently for polymers consisting of composite shells ( Ollila et adl2011al ) like those 



in th is work. The method is similar to Peskin's immersed boundary method (jPeskin 
2002). The LB metho d has been used successfully to model fluid flow in porous me- 
dia (jKang e7aZll2002l ). 

We generate the node distribution for the shell (nsh(r)) in simulation by using the 
atomic coordinates of spherical fullerenes, consisting of N carbons, scaled so the nodes 
sit at radius R. The number of nodes is chosen sufficiently large for a given R to guarantee 
a node placement denser than the resolution of the underlying lattice. This is necessary 
in order to resolve the spherical shape as the value of the coupling constant 7 is increased 
away from the free-draining limit. The uniform-density sphere (|3.1j) and the shell of finite 
thickness (|3.10j) discussed below are generated by placing nodes at intervals Ax n on a 



cubic lattice at coordinates r which fulfill ^ |r — r c 
respectively. 



< R and iJi < r - r cm | < R 2 



3. Analytical Results 

Next, we summarize and present calculations of the components of F and T in typical 
steady states for particles with different node densities n(r) that we will later compare 
to our simulation result. We then look closely at the dynamic problem of the particle 
oscillating sinusoidally in the fluid. 
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3.1. Steady-state solutions 
In this subsection of the paper, we summarize previously derived steady-state solutions 
to (|2.ip . where dtu is assumed to be zero, for spherical particles of different node densities 
that will act as limiting cases of the oscillating particle solution. The particle is assumed 
to be fixed in place and the far- field velocity constant. 

The simplest experimentally relevant case is that of a sphere of uniform node density 
(6 is the Heaviside step function), 



n Sp (r) = XQ(R - r) = N^ttR 3 ) ^(R - r), 



(3.1) 



immersed in a background fluid whose far- field velocity Uoo = u(r — > oo) = Uqb z . As 
mentioned in Section 12.31 N is large enough in the simulations so that the discrete 
nodes are spaced cl oser together than the fl uid mesh spacing to approximate a uniform 
density reasonably. iDebve &: Buechel ( 1948| ) solved the problem, in the context of the 
uniform density sphere being a model for a polymer in solution, for arbitrary values 
of the dimensionless parameter (3 = kR = yf jX/r/R. The product 7A describes the 
strength of coupling between the phases. In simulation, one may lower the node density 
A and increase the coupling parameter 7 while keeping the product the same without 
significantly changing the result s as long as the no de placement is sufficiently dense to 
resolve the shape of the object ( Ollila et al. 2011 bl ). Debye and Bueche presented their 
original solutions as functions of /3 as defined here without any reference to the size of the 
node itself, but only to R and a so-called shi elding length y/rj/JXy), which is also referred 
to as mean pore size by Abade et al. ( 2010al) . Debye and Bueche also related the shieldin g 
length to the slip length in the Navier slip boundary condition ( Debve Sz Buechel 1 1948 ). 
Moreover, Bueche found the drag force on the spherical uniform-density "polymer" to 
be F = Fe z , where 



2/3 2 G (/3) 



6tttiRU 2/3 2 + 3G (/3) ' 
G (/3) = l-(l//3)tanh(/3). 



(3.2) 
(3-3) 



Both the limit of zero, /3 — > 0, and infinite, f3 — > 00, coupling in p.2[) give the intuitive 
results F — > and the Stokes formula F — > Fs = Girrj RUq first derived by Stokes fo r 
an impermeable sphere with a no-slip boundary condition on its surface ( Stokes 1880h ■ 
At low porosity (large /?), (|3.2p can be approximated by 



F 



2/3 2 



2f3 2 



(3.4) 



which differs from (|3.2j) by less than 10% for (3 > 10.9. T his approximation will be 
relevant for what follows so we treat it here in more detail. Sutherland fc Tan] (|l970h 
arrived at (|3.4[) directly by assuming Darcy's law (see text below (|2 - 2|) ) to apply in the 
form 



u(r = Re r ) = gU{)G z = —r/n 2 \7p(r = Re r ) 



(3.5) 



on the surface dB(t). By matching the pressure based on (13.511 an d that based on the 
stream function for a solid sphere (see e.g. Landau fc Lifschita (1987)), they found 1— g = 
2/3 2 /(2/3 2 + 3) and hence they called g the permeation coefficient. This shows that the 
function Gq((3) appears due to the Laplacian term of (|2.16[) included in the Brinkman 
m odel but not in t he Darcy model, which is also apparent from the explicit calculation 
of lFelderhoJ (Il975h . 
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Alternatively, one may place the uniform-density sphere in a flow of constant shear rate 
Q and constrain it not to rotate in which case the particle is subjected to drag torque 
T = Te n , where e„ is a normal vector perpen dicular to the shear plane. Its magnitude 
T = |T| is given by (|Felderhof fc Deutcrj|l975l) 



T 3 3 coth (3 

4nr]R 3 Q ~ + ]P ' 



(3.6) 



We define the Stokes torque (jGoldman et aZ.Ml967h T s = Anr]R 3 Q and find T -> as 
P — > and T — > Ts as j3 — > oo. Equation (|3.6|) is found as the solution to the full 
Brinkman problem from the mean- field theory of Felderhof fc Deutch ( 1975T) or as the 

limit of vanishing hard core for a coated particle ( Cichocki fc Fe lderhof 2009)_. 

Felderhof fc Deutch have written a ser ies of publications (jFelderhof fc Deutchl [l975t 



Deutch fc Felderhof! 1 1975t iFelderhofl 119751 ) on fri ctional properties of dilu te polymer so- 
lutions in which they show how the macroscopic iDebve fc Buechd (|1948j ) results for the 
hydrodynamic frict ion coefficients are obtained as mean-field approximations from a mi- 
croscopic theory bv lKirkwood fc Risemanl (|l948h . The term mean-field is used here as the 
average flow velocity u(r), average pressure p(r) and average force density f (r) are taken 
over the statistical distribution P(ri, . . . , rjy ) of the node positions, w hich were con- 
sidered as segments making up the polymer ( Felderhof fc Deutchlll975h . Felderhof and 
Deutch's work is significant in that they considered more general density distribution 
than nsp(r) in (|3.ip . 



In particular, they considered (jDeutch fc Felderhoflll975t lFelderhoflll975l ) an infinites- 
imally thin shell for which the node density is given by 



nsh(r) = A sh J(r -R) = N^ttR 2 )-^^ - R), 



(3.7) 



where Ash is the uniform surface density. Such shells are of particular interest in biophys- 
ical problems such as leaky vesicles and encapsulated drug delivery. They calculated the 
shell to experience a drag force and torque equal to 



F 2(3 2 
Fs~ ~ 2/3 2 + 9' 
T_ /3 2 
Ts~ ~ /3 2 + 9 



(3.8a) 
(3.86) 



in the same setting that gave (|3.2I) and ([376j) for the uniform-density sphere. We emphasize 
that /3 in (|3.8|) is still equal to Ry^jX/rj, where A from (|3.Ij) guarantees correct units. The 
small difference between (|3.4j) and p.8a| suggests that a low-porosity shell and sphere 
(large (3) might be difficult to distinguish from one another based on drag force. 

We have derived both results of (|3.8p (they are also limiting cases of the dynamic 
calculation we describe in the next subsection) for a stationary particle by requiring the 
force and torque based on the coupling and the fluid stress to match locally at every 
point on the shell, 



-7Ask(-u) = er • e,.; 
-Re r x 7Ash(— u) = Re r x cr ■ e r . 



(3.9a) 
(3.96) 

The fluid velocity field u and stress a that go into (|3.9[) are based on known stream 
functions ( Batchelorlll967 ; Landau fc Lifschita 1987 ) whose constants are left arbitrary 



to be determined by impo sing (|3.9|) . 

Bhatt fc Sachetl (|l994h studied a shell of finite thickness for which the node density 
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can be written as (i?2 > Ri ) 

riAn(r) = A A „ (e(R 2 -r)-Q(Ri- r)) ; 

Rl) 



(3.10) 



We will refer to it as an annulus due to the shape of its cross section. They solved the 
steady state version of (|2.ip in all space by imposing continuity of velocity and shear stress 
at both the inner, r = R%, and outer, r = i? 2 , surface. They found (Fs = 6 7r 77 R 2 Uq, 
Pi = Riy/l^An/v, i = 1.2) 



F/F S 



1 r 



(cosh (3 2 
(sinh /3 2 



sinh P2 
cosh P2 



(3.11) 



where the functions Hi and H2 can be found in lBhatt k, Sacheti (1994). Moreover, the 
right-hand side of (|3.11j) reduces to that of (|3.2j) in the limit /3\ — > by identifying 

3.2. Oscillating particle 

In this section, we study the hydrodynamic force F experienced by the infinitely thin 
shell and the uniform-density sphere oscillating in a fluid along a fixed axis. We briefly 
s ummarize recen t work on the subject and present new time-dependent solutions to (|2.1[) . 



Stokesl (|1901l ) was the first to solve for the hydrodynamic force exerted on a solid 



sphere oscillating sinusoidally at an angular frequency w in a quiesc ent incompre ssible 
fluid. His result and its generalizations are recapitulated by both iLambl (|1932l ) and 
Landau fc Lifschitzl (jl987h . The velocity of the particle is assumed to be v = v e z = 



voe lLUt e z , which, due to symmetry considerations, allows one to describe the resulting 
fluid motion outside the particle as a doublet-stokeslet combination for which the stream 
function ipo reads in spherical coordinates 



4>0 = ho{r) v sin 2 0, r R; 
h (r) =A/r + (D/k)(l + l/(kr))e 



-kr 



(3.12a) 
(3.126) 



where we have used the abbreviations k = (1 + i)a and a = yj pco/2rj, and ho{r) contains 
only terms that vanish as r — > 00. The polar angle 9 is measured with respect to the 
z axis, i.e. direction of particle or ambient flow velocity. Stokes solved the linearized 
equation (|2 . 1 [) by imposing a no-slip boundary condition on the surface of a solid sphere: 
u(|r| = R,6) = v(|r| = R, 9), which corresponds to the limit 7 —> 00 in (|2.5|) . The 
imposition of impenetrability and the no-slip boundary condition on the (outer) surface 
of the sphere amounts to four equations for the complex coefficients A and D in (|3.12|) . 
For any of the particles of (|3.ip , Q3.7P or (|3.10p , the coefficients A and D will in general be 
different due the node distribution inside B{t), but when the coupling constant 7 in (|2.5[) 
goes to zero they should also go to zero for an y value of a. A suitable stre am function 
for solving (|2.1j) inside the sphere has the form ( Vainshtein fc Shapiro 20091 ) 

hi(r) v sin 2 9, r ^ R; 



tpi 
hi(r) 



Br 2 +C 



'sinh(fc/r) 
kir 



cosh(fcjr) 



(3.13a) 
(3.136) 



where kiR= \J fi 2 + 2iY 2 and Y = Ra. If the Laplacian term of (|2.1j) is neglected inside 
the sphere, we may set C to zero in (|3.13|) . We remark that the parameter Y is related 
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to the Womersley number ( Womer sle vl 1 1 9 5 5 ) by W = V2Y, which expresses the ratio 
of oscillatory fluid inertia to the shear force. 

We may proceed to calculate the fluid velocity field components, u r and ug, and com- 
ponents of the stress tensor, o~ r g and o~ rr , in a given region of space via 

u r = (r 2 sin0) -1 d e V; u e = -{r sin 8y 1 d r ip; (3.14) 

fldUr . dug Ug\ 

arf = t¥ + ¥-7) ; (3 - 15) 
o- rr = -p + 2r i-^-- (3-16) 

The pressure, p, is associated with the irrotational part, i.e. ^doublet = (A/r)v sin 2 9 in 
(I3.12flp . of the flow and is solved from Vp = —pdud ou b\et/dt. We note that in determining 
the pressure distribution inside the sphere in the DBB model, the contribution due to 
the irrotational part of the force f must be included. 

Once the stress tensor is known, (|2.3[) or (|2.5|) may be used to calculate the total 
hydrodynamic force on the sphere with arbitrary A and D as 



F = Fe, = -2wpuivR 3 



i 31+ kR 



tle z , (3.17) 



.3 2 Y 2 

where the dimcnsionlcss function O = d Re + i£li m reads 

9{1 + (1 + i)Y) + 2iY 2 ' v ' 

and A and D have units of length cubed and length, respectively. Equation (|3.17[) has 
the attractive feature that d = 1 corresponds to the impermeable no-slip result for a 
sphere. Setting ft = 1 and taking the real part reduces it to (v = voe lult ) 

TZ{F} = -QirriRvo cos(wt)(l + Y) (3.19) 

- (4/3)7rpJ? 3 (-^o)sin(a;t)(- + 

which is the no-slip result valid for a solid sphere obtained bv lStokesl(ll90lh . where TZ{F} 
denotes the real part of F. 

For the more general case (Vl ^ 1), 71{F} can be written as 

K{F} = —QuriRvQ cos{ujt)C & 

~ (4/3)7rpi? 3 (-wv )sin(wi)C A d; (3.20a) 

C s = ( 1 + y) Vi Rc - ( Y + \y 2 ) Im ; (3.206) 



1 . 9\„ f 9 



9 

° Ad = { 2 + 4Y ) nRc + 4^)° Im ' (3 - 2 ° c) 

where we have again used the fact that v = iuiv. The factors C s and Cx& depend both 
on the frequency and the boundary condition on the surface of the sphere (through A 
and D). Having derived (|3.20[) , Vainshtcin and Shapiro attempted to generalize (|3.19|) to 
a uniform-density sphere (sec (|3.1[1 ) by applying the boundary condition (|2.8[) and using 
(|2.7|) inside the sphere without the Laplacian term. They refer to this as the "Darcy 
model" to which their solution reads ( Vainshtein fc ShaDiro|[2009| ) 

2/3 2 + 4»y 2 

VS ~ 2/32 + 3 + 3((1 + i)Y + 2iY 2 ) ' [6 > 
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where B = R^J^X/r}. Equation (|3.21[) is, however, unsatisfactory since 

4iY 2 

Ovs — > t, ; ^rr when B — > 0. 

Vb 3 + 3((l+i)F + 2ir 2 ) H 

When B —¥ the particle is completely permeable and the hydrodynamic force should 
vanish altogether (the fluid and particle phase become completely decoupled in the frame- 
work of the DBB model (|2.1|l ). They also presented an approximate generalization for the 
Brinkman problem with the Laplacian term of (|2.1[) inside B(t) included, but it inher- 
ited the problem of (|3.2ip . We may lift the discrepancy of (|3.2ip by considering the DBB 
model of (|2.1[) without the Laplacian term inside B(t). Thus, we set C=0 in (|3 . 1 3 £»[) and 
use (|3.14[) - (|3.16[) to calculate the resulting velocity and pressure fields inside and outside 
B(t) via p,13|) and (|3.12p , respectively. We require continuity of velocity and pressure 
on the surface r = Re r as opposed to (|2.8[) . This formulation gives 



2B 2 

° D = 2p + 3 + 3((l + i)Y + 2iY*)' (3 ' 22) 

which features the limits £Id — > &s B — > and f2.o — > 1 as — >oo. The reason for the 
unphysical B — > limit of (|3.21j) can thus be attributed to the boundary condition (|2.8j) . 
which does not correspond to continuity of fluid velocity at the boundary if v 0. The 
Y — > limit of (J3T22]) agrees with (|3.4p . This now gives us an expression with physically 
re asonable limits, b u t only within the Darcy model approximation. 



asonable limits, b u t only witnm the Darcy model approximation. 
lLooker fc Carniel ( 20041 ) considered the hydrodynamic force experienced by a sinu- 



soidally oscillating nearly impermeable sphere using a perturbative expansion in 1/B. 
However, they assumed the normal component of the fluid velocity to be zero on the 
surface of the particle, i.e. u • e r = for r E dB(t). They found the hydrodynamic force 
to be 

F/F s = l + k + k 2 /9 - (1 + kf/{BQ, (3.23) 

where B = kR as defined in the present work and £ i s a slip coefficient rela ted to the 
tangential fluid slip length 3 by 3 = l/«) = R/(/3£) (|Looker fc Carniell2004[) . We may 



compare their model to ours by using (|3.20[) to extract the coefficients C s and Ca<i, which 
become 

C s = l + Y _(1 + 2Y)~; (3.24a) 
R 

1 9 9 / 1 \ 5 

C A d = - + — --(! + -)-• (3.246) 

Ad 2 AY 2 V YJ R V ; 

We note that within the DBB model, we do not concern ourselves with a decomposition 
of the hydrodynamic coupling, i.e. 7, to a normal and a tangential component. However, 
in order to compare Looker and Carnie's result to ours, we need to determine a value for 
£. In keeping with their theory, we have assumed £ to be a fixed number, which we have 
determined to be £ = 0.9 by least squares fitting the Y — > limit of (|3.24a[) to (|3.2p in 
the range B G [5, 100]. We will use this value of £ throughout the paper as a similar fit to 
(|3.4|) gives poorer agreement. One immediately sees that the factors in p.24j) reduce to 
those in p,19|) as the slip on the surface vanishes 3 — > 8 — > 00, but due to the nature 
of the expansion, the weak-coupling limit is not correctly reproduced. We comment more 
on (pOl)) in the Results. 

3.3. Exact solution for the oscillating shell 

Finding full non-perturbative solutions to the time-dependent equations (|2.1j) with gen- 
eral force densities like those in (|2.2j) . not just using the Darcy approximation, has not 
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been done. In particular, there does not exist an analytical solution to (j2.1j) for the shell 
of (|3.7p oscillating in a viscous fluid. We establish such a solution here by requiring a 
matching condition between the force on the shell determined both by the coupling 
in the DBB model and the fluid stress 



/ -7"Sh(r)(v - u)d 3 a; = / 

JB(t) JdB(t) 



■AS, 



(3.25) 



where the integrand on the left-hand side is a force per volume. Equation (|3.25p has 
sound physical limits as it guarantees there will be no force acting on the particle if 
7 = since in that case the only solution is to have A = D = 0. Integrating over r 
changes the left-hand side into a surface integral and we may equate the integrands, i. e. 



-7Ash(v - u) = e r • er 

"7^Sh(iV - u r ) — a rr ; 
-j^Sh.{ve —ue) = <J r e 



(3.26) 



is required to hold for ^ 9 ^ tt at r = R at any given time; u and <x are based on the 
stream function of (|3. 12 &[) . By solving (|3.26p . we obtain closed- form expressions for the 
unknowns A = A-r c + iA\ m and D = Z3r c + iDj m in terms of /3 and Y . We note that 
the algebra involved is simplified by the substitution D = D exp kR and by the use of 
a symbolic computation package like Mathematica. Plugging the resulting expressions 
into (|3.18[) gives the real and imaginary parts of (|3.18[) for the shell as fine = and 
f^i m = where 

fii = /3 8 (81 + 162F + 162F 2 + 36F 3 + 4F 4 ) 

+ (9/2)/3 6 ^729 + 1539F + IQ2QY 2 + 486F 3 + 80F 4 + 12F 5 ) 

+ 18^2187 + 5103F + 5832F 2 + 2430F 3 + 504F 4 + 144F 5 + 5AY 6 + 4F 7 ) 

+ 1458^ 2 (1 + Y) (81 + 162F + 162F 2 + 36F 3 + 9Y 4 + 6Y 5 + 2Y 6 ) ; (3.27a) 



fl 2 = -Y ^(3/2)/5 6 (9 + 2Y) [27 + MY + 5AY 2 + 12Y 3 + AY 4 
+ 18^729 + 2F(810 + 891F + 324F 2 + 81F 3 + 21F 4 + AY 5 ) 
+ 162^ 2 (9 + 2Y) (81 + 162F + 162F 2 + 36F 3 + 9F 4 + 6Y 5 + 2Y 6 ' 



(3.276) 



tt 3 = [81 + 162F + 162F 2 + 36F 3 + 4F 4 j (/3 8 + 9/3 6 (5 + Y) 

+ (9/4)^ 4 (297 + 162F + 18F 2 + AY 3 + AY 4 ) 

+ 27^ 2 (135 + 162F + 54F 2 + 12F 3 + 6Y 4 + 2Y 5 ) 

+ 81(81 + 162F + 162F 2 + 36F 3 + 9F 4 + 6Y 5 + 2Y 6 )) , 



(3.27c) 



which is now an exact result for the shell. 

The relevant limits of (|3.27p merit comment. The zero-frequency limit (Y — > 0) for our 
model gives 

F/F s -> 2/3 2 /(2/3 2 + 9) 
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and 

C A d -> 0. 

This is consistent with performing the force matching using the zero-frequency stream 
functions that yields the steady-state shell results of (|3.8[) . The infinite- frequency limit 
gives 

lim F/F s = lim C s = 2/3 2 /9 

and 

lim u>C Ad = 0. 

If we allow j3 to tend to infinity, these four limits agree with those for the no-slip boundary 
condition. For the infinite-coupling limit (/? — > oo, Y given), we find 

c s -> i + r 

and 

C Ad ^ l/2 + 9/(4Y), 

which agree with Stokes's no-slip result of (|3 . 19f) . Thus, f|3 . 2T[) matches up with all the 
known relevant limits. 

3.4. Exact solution to the oscillating uniform- density sphere 

Another case for which there exists no analytical solution to the time-dependent equation 
(|2.ip is that of the uniform-density sphere of (|3.1|) oscillating in the fluid with velocity 
v = i>oe lult e z . We have solved this problem by formulating the system of equations (the 
fields based on (|3.12|l / ()3.13j) are denoted with/without the tilde), 

u r = u r , ue — ug; a rr = cr rr , a r e = ove,V r G dB(t) (3.28) 

which are eight equations for the real and imaginary components of the four complex 
numbers A, B, C and D of (f3~l~26|) and ([3T3&]) . The use of a symbolic computation pack- 
age vastly simplifies the task of deriving these results. The resulting real and imaginary 
parts f^Ro = O1/T23 and Oi m = ^2/^3 are 

fii = 2p 2 (p 1 cos{2Y 2 /Z) - 2ZYP 2 sva{2Y 2 /Z) 

+ P 3 cosh(2Z) + 2ZP 4 sinh(2Z)) ; (3.29a) 



^2 = -&P ( YP 5 cos(2Y z /Z) - 2ZP 6 sm(2Y z /Z) 

+ YP 7 cosh(2Z) + 2ZYP S sinh(2Z)) ; (3.296) 

il 3 = (81 + 162Y + 1Q2Y 2 + 36Y 3 + 4F 4 ) 

({Z 2 - Y 2 )(P 9 cos(2Y 2 /Z) - 2ZYP W sin(2F 2 /Z)) 

+ (Z 2 + Y 2 )(Pu cosh(2Z) + 2ZP l2 sinh(2Z))V (3.29c) 

where \f2Z = \J P 2 + \/ P A + 4^" and the Pj = Pj(/3,Y) are polynomials that arc pro- 
vided in Appendix IB] 

The result (|3. 29[) has the low frequency, Y — > 0, limit 

2^ 2 (1 -p-itaohp) 



F/F s 



2/3 2 + 3(l-/3- 1 tanh i 5) 
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0.001 0.01 0.1 

Ny [10" 6 kg/s] 

Figure 2. Comparison of normalized steady-state drag force F/Fs versus N'y for particles 
with different node density distributions n(r). The lines are the theoretical predictions for the 
shell f (|3.8a|l . solid line), uniform-density sphere f ()3.2p . dashed line) and the annulus ( (|3.1ip . 
dot-dashed line). Taking the thin-annulus limit of (|3.11[1 leads to numerical cancellation errors 
for large N^, but the agreement with simulations is still good. 




Figure 3. Normalized steady-state drag force for a shell multiplied by the radius R at which 
the nodes sit. As the radius is increased to R — 3.3Aa;, the simulation results differ from the 
theory by no more than 0.6% at f3 — 14.4. 

and 

C A d -> 0, 

which agrees with the steady-state result of Debye and Bueche Q3.2p . The (3 — > oo limits 
(impenetrable, no-slip limit) are C s — > 1 + Y and CAd — > 1/2 + 9/(4F), in agreement 
with Stokes's no-slip result of (|3.19[) . 

4. Results: Comparison to simulations 

In this section, we compare theoretical predictions to Lattice-Boltzmann simulations 
for various node densities, for both steady state and as the particle oscillates sinusoidally. 

4.1. Steady-state 

In figure [21 we compare theoretical predictions to LB simulations for the steady-state 
drag force on the shell, (|3.8q|) . uniform-density sphere, ()3.2j) . and on the annulus, p. lip . 
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The results are presented as a function of the product iV7 of the node number and 
coupling parameter as one cannot characterize the annulus in terms of a single f3 that is 
proportional to the radius. The measurement is performed in a cubic simulation box with 
periodic boundary conditions in x and z directions and parallel no-slip walls at y = 
and y — 30 Ax moving at velocity Uoo = Uq e z = (O.OOOlAx/Ai) e 2 . The particle is kept 
immobi le at the centre of the box. We choose this arrangement as it reduces finite-size 
effects (jLiron &: Mochonlll976l ) in the velocity field at least to 1/L 2 from 1/L. 



We find our simulation results to be in excellent agreement with the theoretical predic- 
tions without any fitting parameters. The simulation method is therefore able to differ- 
entiate between a shell, an annulus and a uniform-density sphere without ambiguity. It is 
noteworthy that even a sub-grid thick annulus (R2 — R\ = 0.72Ax) is so closely matched 
with the theoretical prediction. This is most likely due to the fact that the node has a 
compact support on th e fluid mesh which in this case couples a node only to the unit 



cell in which it resides (jAhlrichs fc Dunwedll998t ISmith fc Dennistonll2007t lOllila et 
2011bh . 



Figured] is, however, a more ideal case for the simulations where R/Ax = 3.3 was not 
commensurate with the underlying fluid mesh. In figure [3J we compare shells of different 
radius - node count combinations to the prediction of (|3.8a[) . We observe the simulation 
results to agree with the prediction for all (3 = RyJ^X/r] in the case of the largest 
R = 3.3Ax. The only caveat is the apparent mismatch between theory and simulation 
for the smaller shells at large (3. A change in the surface area per node, 4irR 2 /N, did 
not reduce the mismatch, for which reason we ran the simulation for different radius 
- node count combinations. This suggests that the mismatch is purely a discretization 
effect of immersing an off-lattice spherical shell into a cubic lattice. Figure [3] does not 
just highlight the effect of increasing particle size, but it also indicates the importance 
of incommensurability with the lattice. The two largest particles arc nearly of the same 
size, but for the commensurate particle R/Ax = 3, the lattice effects are emphasized as j3 
increases compared to the incommensurate case R/Ax = 3.3 and the lattice effects again 
decrease for R/Ax = 2.5. This commensurability effect makes it difficult to calculate 
analytically the discretization errors from this sort of immersed boundary simulation. 
However, if we fit the analytical curves to the simulation data by allowing the radius R 
to be a fitting parameter, then we find that the fitted radius differs from the radius at 
which the discrete nodes sit by no more than 0.1 Ax for the cases plotted in figure G2 
thus putting a bound on the discretization errors. 



The drag torque experienced by the uniform-density sphere, (|3.6|) . and by the shell, 
(|3.86p . is more sensitive to changes in the radius R than Stokes drag as the dependence 
is cubic. We have measured the drag torque on spheres and shells of radii 3.3 Ax, 3.0 Ax, 
1.67 Ax and 1.44 Ax. We present the measurements together with theoretical predictions 
in figures Eta) and (b) . Figure HJa) shows the theory without any fitting parameters 
whereas in figure E{b), the radius has been used as the fitting parameter. Even though 
the fitting parameter turns out to be a very small effect, giving a fitted radius 0.1 Ax 
larger than the radius at which the nodes sit, it has a noticeable impact as the torque 
is proportional to R 3 . Moreover, the effect of the compact support of the nodes on the 
fluid mesh has a big effect on small (R w Ax) spheres since its support spreads the nodes 
out to radii between Ax and 2Ax. As with the drag force, the discretization errors for 
the torque can be minimized by u sing an irrational r atio R/Ax, for which the particle is 



incommensurate with the lattice (jOllila et «Lll2011bh 
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Figure 4. Normalized steady-state drag torque on uniform-density spheres (solid symbols) 
and infinitely thin shells (hollow symbols) of different radii and node counts, (a) Measure- 
ment data are shown without any fitting, (b) Normalized measurements are shown with AR in 
Ts = 4ttt](R + AR) 3 Q used as the fitting parameter, which falls between AR — 0.10 Ax and 
0.15 As. 



4.1.1. The limit of impermeability 

We conclude the discussion of the steady-state results by referring to the limit of an 
impermeable particle, i.e. when F/F$ = T/Ts = 1. In the context of simulations, we 
refer to this limit as saturating /3 since it requires 7 to be increased until impermeability 
is reached numerically. One might think that one could simulate a porous particle and 
consider it to be equivalent to an impermeable particle with a smaller effective "hydro- 
dynamic" radius. However, based on (|3.8q|) and (|3.8fej) . one cannot simulate nodes placed 
at radius R for a sub-saturation (3 and claim to be modelling a particle with an effective 
radius R based on, say, Stokes drag F = Q-kijRv as a measurement of the drag torque 
will give a very different value for the effective radius, which can be seen, for instance, 
by writing (|3.86|) as 

T = AnriQR 3 = AmjQ^Rfi 2 ' 3 / '(f3 2 + 9) 1/3 ) 3 . 

Clearly the "effective" hydrodynamic radius from the torque R and from Stokes drag, 
R = 2(3 2 R/ (2/3 2 + 9), are very different numbers unless /3 is large. Therefore, in simulating 
an impenetrable sphere, one must increase j3 until the hydrodynamic radii based on 
different measures all agree to within an acceptable level of tolerance. If one chooses a 
value for /? below saturation, one may only claim to be simulating a porous particle of 
the chosen node density n(r). That is, one should be wary of models using an effective 
hydrodynamic radius as there is more than one way to define such a radius and different 
measures will not generally give the sa me measure. The effective-radius concept has also 
been criticized by Abade et al. ( 2010bJ ). 



In theory, the limit of impermeability requires [3 — > 00. However, it is clear from 
figures [3] and 2] that the drag force and torque are typically slightly larger in simulation 
at large (3 than theory suggests due to discreteness of simulations. Thus, the limit of 
impermeability can be reached in practice in actual simulations for a finite, but still 
quite large, value of (3. Based on figures [3] and 21 particles of different radii and node 
distributions reach the limit at different values of f3. The approach to the limit is dictated 
by the node density and the ratio R/ Ax. A slig htly modified simul ation algorithm may 
also be needed for stability at large values of (3. ( Ollila et al. 2011bl ) 



It is also clear from figure [3] that for a given number of nodes 7Y, a shell will approach 
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1000 



Figure 5. The correction factor C s plotted at Y = aR = 10 for our new shell theory (1 )3. 271 1. 
solid line), the corrected Darcy model f ()3.22)l . heavy dashed line), our new sphere theory ( (13.291) . 
dotted line) and Looker and Carnie's perturbative expansion f (|3.24[) , dot-dashed line), and 
Vainshtein and Shapiro's expression f (|3.21)) , thin dashed line). Our two models new to this work 
and the corrected Darcy result have the correct asymptotic zero and infinite /3 limits. 



the impermeable limit at a lower value of 7 than a uniform density sphere, which might 
make it preferable in some simulations. 

4.2. Oscillating particle 

We first compare the different theoretical predictions for the oscillating sphere case. 
The steady-state uniform-density sphere result for the force, l|3.2[l and Q3.4p . and the 
shell result, ()3.8a)) . are qualitatively very similar. In particular, both the steady-state 
drag force and torque go to zero as /3 goes to zero. One would therefore expect similar 
behaviour from the oscillatory force on both the sphere and the shell. Figure [5] shows 
the correction factor C s to the velocity-dependent part of the hydrodynamic drag force 
(|3.20a[) at Y ~aR—R^J pui / (2rj)~\Q as a function of the dimcnsionlcss coupling parameter 
j3 for Vainshtein and Shapiro's Darcy model p. 21)1 . the corrected Darcy model ()3.22ll . 
Looker and Carnie's nearly impermeable sphere (|3.24|) and our shell p. 27)1 and sphere 
model p. 29)1 . Our theories (solid and dotted lines) and the corrected Darcy model capture 
the correct approach to the limit (3 — > in which the hydrodynamic force must vanish. We 
do emphasize that the existing theories are targeted to model a uniform-density sphere, 
but both the sphere and shell should have similar qualitative behaviour. It is interesting 
that Looker and Carnie's model based on homogenization theory does not lead to the 
peak in C s between /3 = 10 and 100 present in other theories. 

We now compare directly the predictions of the different time-dependent models, after 
substitution into (|3.20)1 . to simulation. In the simulation we set the velocity of the rigid 
test particle to be equal to v(<) = i> cos(wi)e 2 and calculated the force exerted on it. Wc 
placed the test particle at the centre X of the simulation box of size 10R x 10R x 20i? 
and imposed periodic boundary condition in the x and z directions and had p arallel 
no-slip walls in the y direction to reduce finite-size effects (jLiron fc Mochon|[l976l ). The 



particle's velocity can be integrated analytically so we set its position to X(t) = Xo + 
(vo/ui) s'm(ujt))e z . The length scale 1/a is a boundary layer thickness associated with the 
acceleration of the fluid next to the particle. As Y — s- (R/Ax), 1/a approaches Ax and 
viscous effects are confined to a narrow boundary layer on the particle surface. Moreover, 
2 Ax is the smallest length scale in the computational method for which central-difference 
velocity gradients can be computed and thus, we expect deviations between theory and 
simulation due to discreteness of the mesh for Y ^ (R/2Ax). So, to explore large Y, we 
need finer resolution. 

The hydrodynamic force on the particle is characterized by its amplitude max and 
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Figure 6. (a) Phase shift and (b) normalized amplitude of the oscillating hydrodynamic force 
as a function of /3 at Y = 0.187 (or, an angular frequency of u = 27r/(5000Af) = 1.25 MHz). The 
hollow circles correspond to simulations of a shell (7i, JV) = (3.3Aa;=330 nm, 540) with Ax = 100 
nm. The data for a uniform-density sphere (R, 7V)=(3.3A;c, 2247) are plotted as solid circles, (c) 
Phase shift and (d) normalized amplitude of the oscillating hydrodynamic force as a function 
of Y (u was varied to change Y) at fixed ft — 4.0. Simulations in (c) and (d) were performed 
with a finer mesh resolution, Aa; = 77 nm and time step At — 0.59 ns. The hollow squares and 
circles correspond to simulations of shells (R, N)=(4.3Ax, 540) and (9.1Ax, 2252). The data for 
uniform-density spheres (R, N) — (4:.3Ax, 5665) and (9.1Aa;, 61805) are plotted as solid squares 
and circles, respectively. In all plots, lines correspond to our shell theory (solid line), corrected 
Darcy theory (dashed), our theory for the uniform-density sphere (dotted) and Looker and 
Carnie's perturbation theory (dot-dashed) with £ = 0.9. 



phase shift <f> = lim^oo arg(F (t)) — axg(v(t)). These two quantities are easy to evaluate 
since all the linear theories for the force we compare here can be written as 

F(t) = \jo? + b 2 cos(o;t + 0), (4.1) 
a = -(6^ V Rv )C s ; b = (67rr / i?i> )(4/9)r 2 C A d; 
4> = — arccos^a/ \J a? + b 2 ^j , 

for v(t) = vq cos(wi). Since we have available theories for a uniform-density sphere and 
the shell, we have simulated both types of particles. The results for the simulations and 
different theoretical models are shown in figure [U As with the steady state results, we 
see that the simulation results for the uniform-density sphere and the shell differ only at 
moderate to high values of j3 and the amplitude of the force vanishes in both cases as 
0. 

We look first at the phase difference between the set velocity v and the measured hydro- 
dynamic force F as a function of ft = Rk for the angular frequency oj = 2 7r/(5000 At). 
This translates to Y = 0.187 with R = 3.33Ax and 1/a = 17.8Ax for which inertial 
effects in the fluid and compressibility effects at the scale of the particle should be negli- 
gible, but oscillatory effects should still be clearly visible for large j3. We have measured 
the phase shift from simulations and plot the results together with theoretical predictions 
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in figure Eta). The velocity amplitude was set to vq = 10~ 5 v c - In steady state, the drag 
force is opposite to the particle velocity, which means there is a phase shift of — ir radians 
between the two. This is found to carry over to the oscillatory case for small (3. How- 
ever, as j3 increases, more of the fluid around the particle is dragged along with it (i.e. 
its virtual mass increases) and the phase shift drops somewhat. We find all theoretical 
models to give similar results for the phase shift at large /3 corresponding to a nearly 
impermeable particle. However, for a highly porous particle (small /3 corresponding to 
k^ 1 > i?), Looker and Carnie's theory breaks down as expected. Our new model based 
on force matching on the surface agrees with simulation results at all values of /3 for 
both the shell (hollow circles, (R, N)=(3.3Ax, 540)) and a uniform-density sphere (solid 
circles, (R, N)=(3.3Ax, 2247)). The nodes inside the uniform-density sphere were placed 
at intervals of 0.41 Ax. We may therefore conclude that the phase shift is determined 
largely by the drag exerted on the surface for the different types of node distributions. 

Figure Etb) shows the corresponding data for the normalized amplitude of the oscilla- 
tory force. For j3 > 5, our exact solution to the uniform-density sphere agrees with Looker 
and Carnie's model (|3.24[) with £ = 0.9. Our method of fitting their slip coefficient £ thus 
works well for large (3 providing a mapping between the DBB model and the solution 
obtained via the homogenization procedure. Moreover, the Darcy-like model predicts a 
larger force amplitude than the Brinkman-like models, which is consistent with the full 
Brinkman solution having the reduction factor (|3.3|) whose value is between zero and 
unity. The difference between the uniform-density sphere and shell for the simulations 
is also consistent with the results for the steady-state drag force which showed a simi- 
lar difference in figure [2j As /3 approaches zero, the Darcy model of (|3.22j) corrects the 
discrepancy of (|3.21j) that predicts a finite hydrodynamic force on the particle, which is 
unphysical. Our new models of sections 13.31 and [3.41 based on local force matching agree 
well with the simulation data for all f3. As in the steady-state case, the sphere and shell 
results coincide as f3 —> 0. 

Last, we scan over the other parameter, Y = Ra, in the model for a fixed value of 
f3 = R\fynjr\ to see how the phase shift and normalized force amplitude arc affected. We 
could perform simulations for different R such that fixed value Y would be equivalent to 
different values of w. However, this should not make a difference to the force amplitude 
normalized by Fg. Alternatively, if we change Y by adjusting the shear viscosity rj, we 
must also change either 7 or A (or both) to keep (3 unchanged. In addition, a large 
increase in the shear viscosity can make a filled shell effectively a uniform-density sphere 
with a coating. In the end, we fix a few values of R and vary u> to change Y. We pick 
an intermediate value of j3 = 4.0 for which different systems (shell versus spheres) and 
theories are clearly distinguished in figure E^b). A finer resolution, Ax = 77 nm and time 
step At = 0.59 ns was used for these simulations. This was necessary as the length scale 
1/a that appears in ()3.12j) becomes small in lattice units at higher frequencies for the 
original parameters. The length scale 1/a sets the wavelength and decay length of the 
disturbance caused by the moving particle and was too short to be accurately resolved 
on the original mesh. The measurement of the phase turned out to be very sensitive to 
the box size for which reason we used a mesh as large as 12i? x 16R x 28i?. In figure Etc) 
we expectedly find all models to give reasonable phase shifts close to <fi = — ir at low 
frequencies (Y < 0.1). As Y is increased, the phase shifts determined from simulations 
both for the shell (hollow symbols) and the sphere (solid symbols) follow the model 
results well for Y < 1. Again, the phase shift for the shell and sphere are very similar 
and thus appears to be determined on the surface of the particle. 

The normalized force amplitude in figure (6jd) confirms the findings of figure (6jb) 
both by theory and simulation: the uniform-density sphere (filled symbols) experiences a 
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smaller drag force than the shell (hollow symbols) . The simulation data for the uniform- 
density sphere agree well for Y < 0.8 with our theory for the uniform-density sphere 
(dotted line) and with Looker and Carnie's theory with £ = 0.9 (dot-dashed line) (as 
seen for Y = 0.187 and ^ 3 in^b)). The shape of F/F$ as a function of Y for the 
shell, however, is captured quantitatively only by our new shell model for Y < 1.1. The 
simulation data indicate a larger force on the particle than the theory roughly when 
Y ^ 1. Doubling any of the dimensions of the already large simulation volume did not 
change the results either for which reason finite-size effects are not responsible for the 
difference either. We can therefore tentatively suggest that assumptions of the near field 
in linear Stokes theory break down when 1/a ^ R. We therefore restrict the regime 
of validity of the present theories to Y < 1.0. It appears that Looker and Carnie's 
perturbative treatment maps closely to our uniform-density sphere model after fitting 
of £ for Y < 1.0 even for = 4. The effect of fitting becomes less important as is 
increased, a result that is similar to what was stated in IVainshtein fc Shapiro] (|2009h . A 
potentially important practical finding is that for = 4, force matching on the surface 
reasonably captures both the sphere and the shell measurements. The corrected Darcy 
model (dashed line) does have the correct qualitative shape as a function of Y compared 
to the other models and our simulations, but it predicts too large a hydrodynamic force. 



5. Conclusions 

We have validated theoretical predictions of the steady-state drag force and drag 
torque on porous shells of different thickness and uniform-density spheres using Latticc- 
Boltzmann simulations of the full Navier-Stokcs equations. We have found the drag force 
to be in quantitative agreement with the theory without any adjustable fitting parame- 
ters. The torque measurement proved to be more sensitive to discretization effects and 
commcnsurability between the particle and the underlying fluid mesh. 

We derived new closed-form expressions for the hydrodynamic force on an oscillating 
shell and on a uniform-density sphere as a function of the coupling to the fluid. Our 
approach is in good quantitative agreement with simulations and it is consistent for all 
degrees of particle porosity, which is an improvement over existing analytical models. We 
have demonstrated that our models are able to predict the hydrodynamic force correctly 
when the porosity or the frequency of oscillation are changed independently. We have 
also pointed out the regime of validity of these theories that base themselves on the linear 
Stokes equation. 

This work was supported in part by the Academy of Finland through its COMP CoE 
grant, by the Natural Science and Engineering Council of Canada and by SharcNet. We 
also wish to thank CSC, the Finnish IT centre for science, for allocation of computer 
resources. S.T.T.O. would like to thank the National Graduate School in Nanoscicnce, 
the Finnish Cultural Foundation and the Alfred Kordelin Foundation for funding. 



Appendix A. LB method 

Lattice Boltzmann ( LB) is an increasingly common scheme to solve the Navier-Stokes 
equations (jSucci l200ll:ISukop fc ThornelbOOd ) which we summarize here. The method is 
based on solving an approximation of the Boltzmann transport equation (BE) on a cubic 
mesh (or grid) with sites x = (i, j, k)Ax connected to their neighboring sites by a set of 
n vectors {ei}™=Q along which material is transported according to a discretized version 
of the BE. We define a distribution function <?i(x, t) where i labels the lattice directions 
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from site x. For three-dimensional systems, we use a 15-velocity model (Sued 2001) on 
a cubic lattice with lattice vectors e* = (0,0, 0)v c , (±1,0, 0)v c , (0, ±l,0)f c , (0,0, ±l)f c , 
(±1, ±1, ±l)i> c . Physical variables are defined as moments of the distribution functions 
by 

PCM) = ^29i{x,t), (pu a )(x,t) = ^Tgi(x,t)e ia . (Al) 

i i 

The distribution functions evolve in time according to (jBhatnagar et al. 1954 ) 

D i9i = (d t + e ia d a ) gi = --( gi - g- q ) + Wi, (A 2) 



where we have also defined the material derivative Di and a driving term Wi . By choosing 
appropriate moments for the equilibrium distribution g° q and the driving term Wi as 



i 
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(A3) 
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u a Fp + F a up, 



(|2.9p and (|2 . 10[) can be obtained f rom (IA 2[) via a Chapman-Enskog expansion similar 
to derivations in IChen &: Doolen ( 1995 ) . The finite-difference scheme we use to solve 
(|A 2[) is discussed in Ollila et all ( 2011a ) (although thermal noise is not used here). We 
emphasize that the results of the present work are by no means specific to this finite- 
difference algorithm, but it allows a stronger coupling between the fluid and the particle. 
In the present work, the particle is either held fixed or moved sinusoidally in which case 
its equation of motion can be solved analytically and used in the simulation. 



Appendix B. Polynomials Pj in the solution to the uniform-density 
sphere 

The polynomials in (|3.29p are reproduced in full below. 
Pi = (Y 12 + Z 12 )(81 + 162F + 162F 2 + 36Y 3 + 4F 4 )(3 + 3Y + 2/3 2 ) 

+ 2F 4 Z 6 (243 + 729F + 972F 2 + 378F 3 - 324F 4 - 480F 5 - 312F 6 - 216F 7 
+ 2F(27 + 54F + 48F 2 - 36F 4 - 8Y 5 )fi 2 ) 

- F 8 Z 2 (243 + 729F + 972F 2 + 1350F 3 + 1260F 4 + 888F 5 - 96F 6 + 168F 7 

- 2F(81 + 162F + 144F 2 + 12F 3 - 24F 4 - 8Y 5 )(3 2 ) 

- Z 10 (243 + 729F + 972F 2 + 702F 3 + 684F 4 + 1032F 5 + 912F 6 + 264F 7 

- 2F(81 + 162F + 144F 2 - 12F 3 - 48F 4 - 8Y 5 )(3 2 ) 

- 3F 3 Z 8 (36(3 + /3 2 ) - F(-321 - 651F - 732F 2 - 470F 3 + 8Y 4 + UY 5 
+ 2(9 + 50Y + 82Y 2 + 28Y 3 + AY 4 )fi 2 )) 

+ 3F 7 Z 4 (36(3 + P 2 ) + F(159 + 165F + 84F 2 + 74F 3 + 8F 4 + 44F 5 

+ 2(81 + 130F + 98y 2 + 28F 3 + 4F 4 )/3 2 )) (B 1) 
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p 2 = z 10 (27or + 8ioy 2 + losoy 3 + 648r 4 + iosy 5 
+ (8i + isoy + i98y 2 + 72y 3 + sy 4 )^ 2 ) 

+ 4y 4 Z 6 (81 + 153y + 225y 2 + 198y 3 + 204y 4 + 54y 5 

+ (27 + 54y + 6oy 2 + 3oy 3 + 4y 4 )/3 2 ) 

+ y 8 Z 2 (-324 - 234y + 234y 2 + 720y 3 + 264y 4 + losy 5 
+ (27 + 84y + 9oy 2 + 48y 3 + 8y 4 )/3 2 ) 

+ 2y 5 Z 4 (-243 - 729y - 972y 2 - 432y 3 + 204y 4 + 24oy 5 - 36y 6 

- 2y(27 + 54y + 48y 2 + 18y 3 + 4Y 4 )f3 2 ) 

+ y 9 (243 + 729y + 972y 2 + 8ioy 3 + 492y 4 + 324y 5 - 12y 6 

+ (81 + 81Y - 126y 3 - 36y 4 - sy 5 )^ 2 ) 

+ yz 8 (243 + 729y + 972y 2 + 702y 3 + 252y 4 + 6oy 5 - 6oy 6 

+ (135 + 243y + 216y 2 - 18y 3 - 36y 4 - 8y 5 ),3 2 ) (B 2) 



p 3 = (y 12 + z 12 )(8i + i62y + i62y 2 + 36y 3 + 4y 4 )(3 + 3y + 2/3 2 ) 

+ 2y 4 Z 6 (-243 - 729y - 972y 2 - 378y 3 + 324y 4 + 480y 5 + 312y 6 + 216y 7 

+ 2y(-27 - 54y - 48y 2 + 36y 4 + sy 5 )^ 2 ) 

+ Z 10 (243 + 729y + 972y 2 + 1350y 3 + 1260y 4 + 888y 5 - 96y 6 + 168y 7 

+ 2y(-8i - i62y - i44y 2 - i2y 3 + 24y 4 + sy 5 )/3 2 ) 

+ y 8 Z 2 (243 + 729y + 972y 2 + 702y 3 + 684y 4 + 1032y 5 + 912y 6 + 264y 7 

+ 2Y(-81 - 162y - 144y 2 + 12Y 3 + 48Y* + 8Y 5 )f3 2 ) 

+ 3y 7 Z 4 (-36(3 + f3 2 ) + y(-321 - 65iy - 732y 2 - 470y 3 + sy 4 + 44y 5 

+ (is + iooy + i64y 2 + 56y 3 + sy 4 )^ 2 )) 

+ 3y 3 Z 8 (36(3 + I3 2 ) + Y (159 + 165Y + MY 2 + 74Y 3 + 8Y 4 + 44y 5 

+ (162 + 260y + 196y 2 + 56y 3 + 8y 4 )/3 2 )) (B 3) 

Pi = y 11 (27oy + 8ioy 2 + losoy 3 + 648y 4 + iosy 5 
+ (8i + isoy + i98y 2 + 72y 3 + sy 4 )/3 2 ) 

+ 4y 7 Z 4 (81 + 153y + 225y 2 + 198y 3 + 204y 4 + 54y 5 

+ (27 + 54y + 6oy 2 + 3oy 3 + 4y 4 )/3 2 ) 

+ y 3 z 8 (-324 - 234y + 234y 2 + 72oy 3 + 264y 4 + iosy 5 

+ (27 + 84y + 9oy 2 + 48y 3 + 8Y 4 )/3 2 ) 

+ 2y 4 Z 6 (243 + 729y + 972Y 2 + 432Y 3 - 204y 4 - 240Y 5 + 36y 6 

+ (54y + iosy 2 + 96y 3 + 36y 4 + sy 5 )/? 2 ) 

+ y 8 z 2 (-243 - 729y - 972y 2 - 702y 3 - 252y 4 - eoy 5 + eoy 6 

+ (-135 - 243y - 216y 2 + 18Y Z + 36y 4 + 8Y 5 )f3 2 ) 

+ Z 10 (-243 - 729y - 972y 2 - 810y 3 - 492Y 4 - 324Y 5 + 12y 6 

+ (-81 - 81Y + 126y 3 + 36y 4 + 8Y 5 )/3 2 ) (B 4) 
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P 5 = (Y 12 + Z 12 )(81 + 324F + 486Y 2 + 360y 3 + 76Y 4 + 8Y 5 ) 

+ 2F 4 Z 6 (81 + 180F + 270y 2 + 216Y 3 + 136Y 4 - 72Y 6 - 16Y 7 

+ 2(9 + 18Y + 20y 2 + 16Y 3 + 12Y 4 )(3 2 ) 

+ Y 8 Z 2 (-81 - 72Y + 270 Y 2 + 648Y 3 + 536F 4 + 104F 5 - 40F 6 - 16Y 7 
+ 2(27 + 36Y + 2AY 2 - 12Y 3 + 8Y 4 )(3 2 ) 

- Z 1Q (81 + 288F + 378y 2 + 216Y 3 + 401' 4 + 104F 5 + 104F 6 + 16Y 7 

- 2(27 + 72y + 96F 2 + 60F 3 + 16Y 4 )(3 2 ) 

+ YZ S (108 + 32AY + 432F 2 + 63F 3 - 108F 4 + 58Y 5 + 440F 6 + 164F 7 + 24F 8 
+ AY (9 + 18F + 16y 2 - 4y 3 - AY 4 )/3 2 ) 

+ y 5 Z 4 (-108 - 324y - 432y 2 - 225y 3 + 324y 4 + 698y 5 + 568y 6 + 164y 7 + 24y 8 

- 4y(9 + i8y + i6y 2 + 4y 3 + 4y 4 )/? 2 ) (B5) 



P 6 = Y 8 Z 2 (108 + 324y + 432y 2 + 306y 3 + 144y 4 + 76y 5 + 44y 6 + sy 7 

- (45 + 99y + losy 2 + 38y 3 + i2y 4 )^ 2 ) 

+ y 3 z 8 (8i + 288y + 378y 2 + 2iey 3 - 32y 4 - 36y 5 - sy 6 
-(9 + 8y-2y 2 -8y 3 )^ 2 ) 

+ z w (9oy 3 + 2ooy 4 + 22oy 5 + 76y 6 + sy 7 - 3(9 + 27y + 36y 2 + 22y 3 + ay 4 )p 2 ) 
+ y 11 (8i + 252y + 27oy 2 + 72y 3 - io4y 4 - 36y 5 - sy 6 - 3(9 + iey + i4y 2 )^ 2 ) 

- 2y 7 z 4 (8i + i62y + 2i6y 2 + i44y 3 + iooy 4 + 36y 5 + sy 6 

+ (18 + 36y + 28y 2 - AY 3 )(3 2 ) 

_ 4y 4 Z 6 (27 + 81Y + 108Y 2 + 45y 3 - 34y 4 - 58y 5 - 30Y 6 - 4y 7 

+ y(9 + i8y + 22y 2 + ey 3 )^ 2 ) (B6) 



P 7 = (y 12 + Z 12 )(81 + 324y + 486y 2 + 360Y 3 + 76F 4 + 8Y 5 ) 

- 2Y 4 Z 6 (81 + 180y + 270 Y 2 + 216Y 3 + 136Y 4 - 72Y 6 - 16Y 7 
+ 2(9 + 18F + 20Y 2 + 16Y 3 + 12Y 4 )^ 2 ) 

+ Z 10 (81 + 72Y - 270Y 2 - 648F 3 - 536F 4 - 104F 5 + 40F 6 + 16F 7 

- 2(27 + 36F + 24y 2 - 12y 3 + 8F 4 )/3 2 ) 

+ F 8 Z 2 (81 + 288F + 378F 2 + 216F 3 + 40F 4 + 104F 5 + 104F 6 + 16F 7 

- 2(27 + 72Y + 96y 2 + 60F 3 + 16F 4 )^ 2 ) 

+ F 5 Z 4 (108 + 324F + 432F 2 + 63F 3 - 108F 4 + 58F 5 + 440F 6 + 164F 7 + 24F 8 
+ 2Y(9 + 18Y + 16Y 2 - AY 3 - 4F 4 )^ 2 ) 

- FZ 8 (108 + 324F + 432y 2 + 225y 3 - 324y 4 - 698F 5 - 568F 6 - 164F 7 - 24F 8 

+ 4r(9 + i8r + i6r 2 + 4r 3 + 4y 4 )/3 2 ) (B7) 
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P 8 = YZ 8 (108 + 324Y + 432Y 2 + 306Y 3 + 144Y 4 + 76Y 5 + 44Y 6 + 8Y 7 

- (45 + 99Y + 108Y 2 + 38Y 3 + 12Y 4 ),3 2 ) 

- Y S Z 2 (81 + 288Y + 378Y 2 + 216Y 3 - 32Y 4 - 36Y 5 - 8Y 6 

- (9 + 8y - 2y 2 - sy 3 )/? 2 ) 

+ y 9 (9oy 3 + 2ooy 4 + 22oy 5 + 76Y 6 + sy 7 

- (27 + 81Y + 10SY 2 + 66y 3 + 12Y 4 ),3 2 ) 

- ^0(81 + 252y + 27oy 2 + 72y 3 - i04y 4 - 36y 5 - sy 6 - 3(9 + i6y + uy 2 )/? 2 ) 
+ 2y 4 z 6 (8i + i62y + 2i6y 2 + i44y 3 + iooy 4 + 36y 5 + sy 6 

+ (is + 36y + 28y 2 - 4y 3 )/? 2 ) 

- 4y 5 z 4 (27 + 8iy + losy 2 + 45Y 3 - 34y 4 - 58Y 5 - 3oy 6 - 4y 7 

+ (9y + i8y 2 + 22y 3 + 6y 4 ),3 2 ) (B8) 

P 9 = (Z 1Q - y 10 )(9 + 1SY + 18Y 2 + 36Y 3 + 36Y 4 + I2f3 2 + 12Y/3 2 + 4/3 4 ) 

- z s (9 + i8y + 9y 2 + i8y 3 + 54y 4 + 36y 5 + 36y 6 

- i2y(i + 3y + 3y 2 )/3 2 + ay 2 p 4 ) 

+ Y 6 Z 2 (9 + 18Y + 9Y 2 + 18Y 3 - 90Y 4 - 108Y 5 + 36Y 6 

- 12Y(1 + 3y - Y 2 )0 2 + 4Y 2 ^ 4 ) 

+ y 4 z 4 (9 + i8y + i8y 2 - 36y 3 - 2i6y 4 - 72y 5 - 72y 6 

- 12Y(3 + 2Y - 2Y 2 )/3 2 - 8Y 2 /3 4 ) 

- Y 2 Z 6 (9 + 18Y + 18Y 2 + 108Y 3 + 216Y 4 + 72Y 5 - 72Y 6 

+ 12Y(1 - 2Y - 2Y 2 )(3 2 - 8Y 2 (3 4 ) (B 9) 

P w = 2Z 8 (9Y + 18Y 2 + 18Y 3 + 18Y 4 + 2^ 4 ) 

- 3Y 5 Z 2 (3 + 6Y + 24Y 2 + 24Y 3 - 12Y 4 + 2(1 - 4Y - 2Y 2 )/3 2 ) 
+ 3YZ 6 (3 + 6Y + 12Y 2 + 24Y 3 + 12Y 4 + 2(3 + 4Y + 2Y 2 )/? 2 ) 

- Y 7 (9 + 18Y + 18Y 2 + 36Y 3 - 36Y 5 + 6(1 - 2Y 2 )^ 2 - 4Y/3 4 ) 

+ Y 3 Z 4 (9 + 18Y - 36Y 2 + 36Y 4 + 72Y 5 + 6(3 + 8Y + 2Y 2 )0 2 + 8Y/3 4 ) (B 10) 

Pn = (Y 10 + Z w ){9 + 18Y + 18Y 2 + 36Y 3 + 36Y 4 + 12(1 + Y)/3 2 + 4/3 4 ) 
+ Y 6 Z 2 (9 + 18 Y + 9Y 2 + 18 Y 3 + 54Y 4 + 36Y 5 + 36 Y 6 

- 12Y(1 + 3Y + 3Y 2 )/3 2 + 4Y 2 /? 4 ) 

+ Z 8 (9 + 18Y + 9Y 2 + 18Y 3 - 90 Y 4 - 108Y 5 + 36Y 6 

- 12Y(1 + 3Y - Y 2 )P 2 + 4Y 2 /3 4 ) 

- Y 2 Z 6 (9 + 18Y + 18Y 2 - 36Y 3 - 216Y 4 - 72Y 5 - 72Y 6 

- 4Y(9 + 6Y - 6Y)/3 2 - 8Y 2 ^ 4 ) 

- Y 4 Z 4 (9 + 18Y + 18Y 2 + 108Y 3 + 216Y 4 + 72Y 5 - 72Y 6 

+ 12Y(1 - 2Y - 2Y 2 )(3 2 - 8Y 2 /? 4 ) (Bll) 
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P12 = 2F 9 (9F + 18F 2 + 18Y 3 + 18Y 4 + 2/3 4 ) 

+ 3Y 2 Z 6 (3 + 6Y + 2AY 2 + 2AY 3 - 12F 4 + 2(1 - AY - 2Y 2 )0 2 ) 

- 3Y 6 Z 2 (3 + 6Y + 12Y 2 + 2AY 3 + 12F 4 + 2(3 + AY + 2Y 2 )0 2 ) 

- Z 8 (9 + 18Y + 18Y 2 + 36F 3 - 36F 5 + 6(1 - 2Y 2 )/3 2 - AY 

+ F 4 Z 4 (9 + 18Y - 36F 2 + 36F 4 + 72Y 5 + 6(3 + 8Y + 2Y 2 )f3 2 + 8F/3 4 ) (B 12) 
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